function VAR = doProxySVAR_single_EV(VAR,DATASET)
% Estimate VAR on full sample.  
% Identify shocks using subsample with non-missing instrument

X      = VAR.X;              % Y = X*beta + U
Xr     = VAR.X(VAR.I_m,:);   % restricted subsample with non-missing instrument
Y      = VAR.Y;
Yr     = VAR.Y(VAR.I_m,:);
VAR.k  = 1;
VAR.mr = VAR.m(VAR.I_m,:); % m_tilde
[VAR.T,VAR.n] = size(Y);
[VAR.Tr,VAR.nr] = size(Yr);

% Run VAR on longest possible sample
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
VAR.bet   = [X ones(length(X),1)]\Y; 
VAR.res   = Y-[X ones(length(X),1)]*VAR.bet;
VAR.Sigma = (VAR.res'*VAR.res)/(VAR.T-VAR.n*VAR.p-VAR.N);     % -30 for FE?
VAR.resr  = Yr-[Xr ones(length(Xr),1)]*VAR.bet;
VAR.Sigmar= (VAR.resr'*VAR.resr)/(VAR.Tr-VAR.nr*VAR.p-VAR.N); % -VAR.N=30 for FE?

% Instrument = residual regression on X
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z         = VAR.mr;

% Identification with ext instrument m on RESTRICTED sample
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Phib    = [ones(length(Z),1) Z]\VAR.resr; % Regress RF residuals on instrument and constant  
Phib    = Phib(2:end,:);                            % drop constant
VAR.Phib=Phib;
Phib11  = Phib(1:VAR.k,1:VAR.k);                    % First stage : HHD residual on instrument
Phib21  = Phib(1:VAR.k,VAR.k+1:VAR.n);              % Reduced form: lnY residual on instrument
b21ib11 = (Phib11\Phib21)';                         % RF/FS = s(21) / s(11)  
Sig11   = VAR.Sigmar(1:VAR.k,1:VAR.k);
Sig21   = VAR.Sigmar(VAR.k+1:VAR.n,1:VAR.k);
Sig22   = VAR.Sigmar(VAR.k+1:VAR.n,VAR.k+1:VAR.n);
ZZp     = b21ib11*Sig11*b21ib11'-(Sig21*b21ib11'+b21ib11*Sig21')+Sig22;    % Q in Gertler Karadi footnote 4
b12b12p = (Sig21- b21ib11*Sig11)'*(ZZp\(Sig21- b21ib11*Sig11));            % s(12)s(12)' in Gertler Karadi
b11b11p = Sig11-b12b12p;                                                   % s(11)^2.  Identifies s(11)=s^p up to a sign restriction
b11     = sqrt(b11b11p);                                                   % s(11): impact of unit structural HHD shock on HHD/GDP(t)
VAR.b1  = [b11; b21ib11*b11];                                              % s(21): impact of unit structural HHD shock on lnY(t)                                   
VAR.b21ib11=b21ib11;  


%Reliability
%  Sigmm = VAR.m'*VAR.m/VAR.T;
ED    = eye(VAR.k)*sum(sum(Z,2)~=0)/VAR.Tr;
mu1   = Z'*VAR.resr(:,1:VAR.k)/VAR.Tr;
 % PhiPhip = mu1*inv(b11b11p)*mu1';
 % VAR.RM    = inv(Sigmm)*PhiPhip*inv(ED)
Bi = [1/b11+(Sig21-Sig11*b21ib11)'*inv(ZZp)/b11*b21ib11 -(Sig21-Sig11*b21ib11)'*inv(ZZp)/b11];
VAR.et = (Bi*VAR.resr')';

PHI = mu1/b11;
GAM = inv(ED)*PHI;
E  = GAM*VAR.et(sum(Z,2)~=0);
V  = Z(sum(Z,2)~=0)-E;
VAR.RM = inv(E'*E+V'*V)*E'*E;


% Impulse Responses
%%%%%%%%%%%%%%%%%%%%
irs(VAR.p+1,:) = VAR.b1(:,1)/VAR.b1(1,1);
 for jj=2:VAR.irhor
 lvars = (irs(VAR.p+jj-1:-1:jj,:))';
 irs(VAR.p+jj,:) = lvars(:)'*VAR.bet(1:VAR.p*VAR.n,:);     
 end
VAR.irs = irs(VAR.p+1:end,:); 


end



